Nonlinear  Inversion  from  Nonlinear  Filters  for  Ocean  Acoustics 


Robert  I.  Odom 
Applied  Physics  Laboratory 
College  of  Ocean  and  Fisheries  Sciences 
University  of  Washington 
1013  N.E.  40th  Street 

phone:  (206)685-3788  fax:  (206)543-6785  email:  odom@apl.washmgton.edu 

Award  #  N00014-05-1-0087 


LONG  TERM  GOALS 

The  long  term  goals  of  this  research  are  to  develop  practical  and  efficient  algorithms  for  application  the 
nonlinear  inversion  problems  encountered  in  ocean  acoustics.  Such  algorithms  would  be  used  for 
estimating  or  accounting  for  the  effects  of  the  environment  on  acoustic  propagation,  detection  and 
tracking  in  shallow  water. 

OBJECTIVES 

The  specific  objectives  of  this  research  are  adapt  a  specific  nonlinear  filter,  known  as  a  Daum  filter,  for 
acoustic  inversion  of  shallow  water  environmental  properties,  and  to  assess  the  performance  of  this 
nonlinear  filter  relative  to  local  linear  inversion  on  the  one  hand  and  global  methods,  e.g.  Monte  Carl 
methods  on  the  other  hand. 

APPROACH 

Many  inverse  problems  of  interest  in  ocean  acoustics  are  intrinsically  nonlinear,  e.g.  inverting 
measured  pressure  data  for  bottom  and  scattering  properties.  The  solution  to  the  nonlinear  inversion 
problem  is  usually  approached  in  one  of  two  ways.  The  first  way  is  to  assume  a  starting  model,  which 
one  hopes  is  near  to  the  true  model,  then  recursively  solve  a  linearized  version  of  the  inverse  problem 
for  corrections  to  the  starting  model  and  model  covariance.  The  advantage  of  this  approach  is  that  the 
numerical  implementation  of  the  solution  algorithm  is  relatively  straightforward  and  in  a  linear 
problem  the  statistical  properties  are  well  defined  and  will  remain  gaussian  if  they  start  out  gaussian. 
However  linearization  of  a  nonlinear  system  can  produce  biased  estimates  for  two  reasons:  1. 
Linearization  of  the  system  and/or  measurement  equations  may  not  be  a  good  approximation,  and  2. 
Nonlinear  systems  do  not  maintain  gaussian  statistics  as  they  evolve  even  if  they  are  initially  gaussian. 
Another  problem  with  linearizing  a  nonlinear  system  is  that  with  a  poor  starting  guess  the  solution 
algorithm  may  never  converge  to  the  true  answer.  If  the  starting  model  represents  a  point  near  a  local 
minimum  of  the  solution  space,  the  final  solution  will  be  trapped  in  that  local  minimum,  and  never 
converge  to  the  true  answer.  This  can  be  circumvented  by  using  Monte  Carlo  techniques  to  randomly 
sample  the  solution  space  for  starting  models. 

The  other  class  of  solution  methods  attack  the  nonlinear  problem  directly  by  using  simulated 
annealing  or  genetic  algorithms.  The  disadvantage  of  these  directly  nonlinear  methods,  is  that  there  is 
no  way  to  conveniently  propagate  the  statistical  properties  of  the  solution  through  to  the  final  result. 
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One  solution  to  this  problem  is  to  find  the  global  minimum  in  the  solution  spaee,  if  one  exists,  then 
linearize  about  the  solution  representing  the  global  minimum  and  do  a  statistieal  analysis  about  that 
solution.  This  was  done  by  Potty  et  al.(2000),  who  employed  a  genetie  algorithm  followed  by  linear 
analysis  about  the  solution  determined  by  the  genetie  algorithm. 

The  reeursive  algorithms  eommonly  employed  for  the  estimation  of  the  model  and  eovarianee  relative 
to  some  initial  starting  values  bear  a  strong  resemblanee  to  Kalman  filters,  whieh  are  eommonly 
employed  in  target  traeking  algorithms.  The  original  Kalman  filter  was  derived  for  strietly  linear 
systems.  However,  the  Extended  Kalman  Filter  ean  be  applied  to  systems  whieh  are  weakly  nonlinear. 
In  the  late  1980s  Frederiek  Daum,  a  mathematieian  working  at  Raytheon  Corporation,  developed  a 
fully  nonlinear  formulation  to  the  filtering  problem  for  target  traeking  (Daum,  1985,  1986,  1987).  His 
theory  is  elegant,  but  impraetieal  from  an  implementation  point  of  view.  Sometime  later  Sehmidt 
(Sehmidt,  1993)  sueeeeded  in  deriving  an  approximate  algorithm  based  on  Daum's  original  theory, 
and  developed  a  sueeessful  numerieal  implementation  of  a  nonlinear  filter  that  was  a  signifieant 
improvement  to  the  Kalman  and  Extended  Kalman  filters  for  the  type  of  traeking  problem  Sehmidt 
was  interested  in. 

This  researeh  aims  to  develop  an  oeean  aeoustie  inversion  algorithm  based  on  Sehmidfs  (1993) 
implementation  of  Daum's  nonlinear  filtering  theory.  The  purpose  is  to  be  able  to  earry  along  the 
statisties  of  the  geoaeoustie  model  parameters  through  the  inversion  proeess.  The  work  is  mueh  more 
than  a  straightforward  "rename  the  variables  and  eode  it  up".  However,  the  traeking  algorithms  do 
bear  resemblanee  to  the  iterative  inversion  algorithms  for  updating  model  parameter  means  and 
eovarianees  of  an  iterated  inverse  problem  as  in,  for  example,  Menke  (1983).  Most  estimation 
problems  ean  be  east  into  an  interative  form,  whereby  the  state  veetor,  whieh  in  our  ease  is  the  oeean 
bottom  model  veetor,  is  updated  sequentially  as  data  is  added.  Estimation  filter  formulations  are  also 
natural  for  range  dependent  or  time  dependent  environments.  Daum's  original  theory  and  Sehmidfs 
praetieal  implementation  assumes  nonlinear  dynamies  and  a  linear  relationship  between  the 
measurement  and  state  veetors.  In  our  ease  the  measurement  veetor,  eomplex  pressure  say,  and  the 
state  veetor,  the  bottom  model,  are  not  linearly  related.  The  filter  needs  to  be  re-derived  from  serateh 
with  the  measurement  to  state  veetor  relationship  appropriate  for  our  oeean  aeoustie  applieation.  Onee 
re-derived,  it  will  need  to  be  eoded,  and  eheeked  against  results  for  linear  inverse  problems.  Dosso 
(e.g.  Dosso  and  Wilmut,  2002)  at  the  University  of  Vietoria  has  developed  a  Monte  Carlo  inversion 
method  for  the  oeean  aeousties  problem.  This  is  eomputationally  very  intensive,  but  he  does  get  the 
full  probability  density  funetion  (pdf)  for  the  model  parameters.  Beeause  the  Sehmidt  implementation 
of  the  Daum  theory  propagates  the  additional  terms  in  the  mean  and  eovarianee  of  the  state  veetor  pdf, 
it  falls  in  between  the  standard  linear  inversion  methods  and  Dosso's  Bayesian  Monte  Carlo  methods. 

Currently  employed  algorithms  for  nonlinear  problems  sueh  as  simulated  annealing  and  genetie 
algorithms  have  no  meehanism  for  propagating  the  statisties.  What  the  nonlinear  filter  algorithm  will 
do  is  provide  a  natural  meehanism  for  updating  the  statisties  as  a  solution  is  determined.  A  eomparison 
of  the  filter  with  an  algorithm  sueh  as  simulated  annealing  would  be  illuminating,  and  a  valuable 
eheek  on  the  filter  algorithm  itself. 

Filter  type  algorithms  are  ideally  suited  to  inverse  problems  with  time  dependent  oeeanography  or 
range  dependenee.  We  do  not  antieipate  attempting  to  inelude  time  dependent  oeeanography  at  this 
time,  but  we  would  like  to  look  at  the  issue  of  range  dependent  inversion.  The  idea  would  be  to 
sequentially  update  parameter  estimates  as  a  funetion  of  range.  Also  note  that  any  inversion  algorithm 
ean  be  east  into  a  filter  like  algorithm  by  supplying  the  data  sequentially  and  updating  the  model 


parameter  estimates  sequentially  as  data  is  added  to  the  problem,  or  a  smoother  by  eonsidering  the 
eomplete  data  set,  and  working  both  forwards  and  backwards  through  the  data  set.  In  the  end, 
probably  the  best  formulation  to  use  for  a  given  inverse  problem  depends  on  the  noise  statistics.  This 
is  also  something  we  propose  to  investigate. 

Linear  inverse  problems  admit  the  construction  of  both  data  and  model  resolution  matrices.  These 
resolution  matrices  can  be  used  as  metrics  with  which  to  estimate  model  uniqueness  and  data 
predictability.  We  will  be  able  to  construct  resolution  matrices  for  the  nonlinear  problem  and  compare 
them  with  their  fully  linear  equivalents. 

WORK  COMPLETED 

As  stated  above  an  inverse  problem  can  be  recast  as  a  filtering  problem.  A  strictly  linear  problem 
becomes  a  Kalman  Filter  (KF).  A  problem  that  has  been  locally  linearized  becomes  an  Extended 
Kalman  Filter  (EKF),  and  the  fully  nonlinear  problem  can  be  represented  as  a  Daum  -  Schmidt  Eilter. 
We  have  completed  working  computer  code  for  a  simple  EKE,  an  EKE  smoother,  an  iterated  Kalman 
filter  and  an  interated  Kalman  smoother  with  which  to  compare  more  complicated  inverse  models. 

Recent  papers  in  the  geoacoustics  literature  (Dosso  and  Nielson  2002)  show  estimated  probability 
density  functions  (PDEs)  of  ocean  geoacoustic  parameters  via  a  Monte  Carlo  method,  and  the  PDEs  for 
some  parameters  are  far  from  Gaussian  -  some  are  bimodal,  some  are  greatly  skewed.  If  linearization 
methods  are  used  on  such  problems  then  the  resulting  maximum  likelihood  estimate  and/or  variance 
about  it  may  be  misleading  descriptors  of  the  true  solution.  While  the  Monte  Carlo  methods  offer  a 
reliable  way  to  address  such  non-Gaussian  statistics  of  the  inverse  solution,  they  are  very  slow  and 
may  not  convey  an  intuitive  understanding  of  these  statistics  that  an  analytic  expression  might.  Our 
work  this  year  has  also  been  to  continue  researching  analytic  means  to  address  these  non-Gaussian 
statistics  and  their  effects  on  resolution  in  nonlinear  inverse  problems  such  as  the  ocean  geoacoustic 
problem. 

RESULTS 

Recent  work  in  the  filter  community  (Daum,  1986)  aims  to  analytically  map  higher  order  moments  of 
the  state  variable  PDEs  for  the  solution  of  a  nonlinear  problem,  given  a  particular  form  of  PDE.  There 
are  close  connections  between  filter  estimation  theory  and  geophysical  inverse  theory,  and  we  have 
been  working  to  adapt  the  filter  theory  moments  work  to  problems  in  ocean  geoacoustic  estimation. 
This  year  in  the  ONR  nonlinear  inversion  project,  the  relation  between  filter  theory  and  inverse  theory 
for  these  problems  has  been  explored.  Eigure  1  shows  how  the  nonlinear  smoother’s  solution  steps 
through  the  space  of  candidate  solutions,  as  it  homes  in  on  the  objective  surface  minimum  at  the 
solution  point.  The  geophysicists’  iterated  Gauss-Newton  method  takes  this  same  path.  In  the  simple 
synthetic  example  problem  shown  in  Eigure  1,  noisy  wave  arrival  times  “recorded”  at  20  receivers  (the 
white  triangles)  are  used  to  estimate  the  source  location  (yellow  diamond)  in  a  constant  medium;  the 
contours  are  the  objective  function  corresponding  to  the  data  misfit.  Eigures  2a  and  2b  show  error 
plots  (estimated  solution  minus  “true”  synthetic  model)  for  two  extended  Kalman  filter  (EKE)  and 
iterated  Kalman  smoother  (IKS)  runs.  The  run  in  2a  and  the  run  in  2b  correspond  to  two  different 
orderings  of  the  same  receiver  data,  including  the  same  instantiation  of  the  random  noise,  so  the  data 
point  order  is  the  only  difference  between  the  runs.  The  horizontal  axes  in  these  plots  is  the  number  of 
receivers  from  which  data  was  brought  in,  with  an  implied  ordering.  Note  between  2a  and  2b  that 
when  the  ordering  of  the  receiver  data  brought  in  changes,  the  EKE  result  changes  but  not  the  smoother 


result.  The  smoother’s  result  is  a  flat  line  for  receiver  location  because  unlike  in  the  EKF,  the  source 
location  is  made  to  be  constant  regardless  of  the  order  that  the  data  were  read  in. 
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Figure  1:  Path  of  nonlinear  inversion  iteration  steps  through  the  objective  space. 
From  the  initial  solution  estimate  (and  its  associated  uncertainty)  at  the  green  dot 
the  Gauss-Newton  process,  and  the  smoother,  take  steps  (the  dotted  line)  toward 
source  location  (yellow  diamond).  The  triangles  are  receiver  locations.  The 
contours  reflect  the  problem ’s  objective  surface,  which  is  not  paraboloidal  because 
this  2D  source  location  problem  is  nonlinear. 


The  filter  literature  used  thus  far  assumes  a  special  form  of  PDF  for  the  nonlinear  problem,  but 
ultimately  one  would  like  to  generalize  this  PDF  form  much  further  to  cover  a  wide  variety  of  ocean 
geoacoustic  problems.  This  year  a  Monte  Carlo  simulation  is  being  developed  to  explore  the  relation 
between  this  PDF  form  and  the  “true”  statistics  for  an  ocean  geoacoustic  problem  based  on  the 
DFMUS  experiment  at  the  Malta  Plateau.  Uncertainty  and  resolution  aspects  of  this  DFMUS 
experiment  problem  have  also  been  initially  explored  for  multiple  receivers  in  summertime  work  on 
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Figures  2a  &  2b:  Error  plots  (estimated  solution  minus  “true”  synthetic  model)  for 
two  extended  Kalman  filter  (EKF)  and  iterated  Kalman  smoother  (IKS)  runs  on  the 
same  problem  shown  in  Figure  1.  The  run  in  (a.)  and  the  run  in  (b.)  correspond  to 
two  different  orderings  of  the  same  receiver  data,  including  the  same  instantiation 
of  the  random  noise,  so  the  data  point  order  is  the  only  difference  between  the  runs. 
The  IKS’s  result  is  a  flat  line  for  receiver  location  because  unlike  in  the  EKF,  the 
source  location  is  made  to  be  constant  regardless  of  the  order  that  the  data  were 
read  in.  In  clarification,  the  iteration  steps  in  Figure  I  are  different  than  the  steps 
along  the  horizontal  axis  here.  Both  the  EKF  and  IKS  curves  were  iteratively 
updated,  and  the  IKS  results  shown  here  are  associated  with  the  final  solution  point 

in  Figure  I. 


the  ONR  ARL  project  managed  by  John  Tague,  and  presented  at  the  FUSION06  conference  in 
Florence,  Italy  (Pitton,  Ganse,  Krout,  Anderson,  2006).  In  this  work,  a  synthetic  demonstration 
inversion  of  ocean  bottom  properties  based  loosely  on  the  geometry  of  the  DEMUS  experiment 
quantified  the  improvement  via  decreased  uncertainty  and  increased  resolution  of  bottom  loss  and 
scattering  functions  given  the  addition  of  a  second  receiver.  See  Figures  3  and  4  for  part  of  the 
demonstration.  The  Monte  Carlo  part  of  this  work  is  a  little  bit  different  formulation  than  the  synthetic 
demonstration  problem  -  the  demonstration  problem  for  the  FUSION06  paper  estimated  the  full 
continuous  bottom  loss  and  scattering  functions,  which  are  parameterized  by  many  parameters,  thus 
making  Monte  Carlo  simulation  computationally  intractable  in  the  scope  of  this  work.  So  for  the 
Monte  Carlo  simulation  the  problem  is  simplified  to  a  small  number  of  parameters  representing  the 
properties  of  the  ocean  bottom.  This  is  a  strong  regularization  of  the  problem,  based  on  additional 


prior  information  in  the  form  of  cores  and  samples  from  the  area  of  interest.  Finally,  in  future  work, 
the  filter  based  developments  will  be  applied  to  this  same  parameter-estimation  problem  to  compare  to 
the  Monte  Carlo  results. 


Figure  3:  Standard  deviations  of  inverted  bottom  loss  functions  of  grazing  angle  in  the  FUSION06 
demonstration  problem,  showing  the  improvements  via  decreased  uncertainty  when  adding  another 
receiver  into  the  inverse  problem.  Standard  deviation  curves  over  angle  are  compared  for  each 
single-receiver  inversion  and  the  combined-receiver  inversion  case.  Note  that  the  benefits  seen  in 
the  individual  receiver  results  are  essentially  combined  when  both  receivers  are  used  -  each  receiver 
provides  for  a  low  standard  deviation  in  a  particular  angular  range,  which  together  provides  for  low 

standard  deviations  over  a  much  wider  angular  range. 
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Figure  4:  Resolution  of  inverted  bottom  loss  functions  of  grazing  angle  in  the  FUSION06 
demonstration  problem,  showing  the  improvements  via  increased  resolution  (equivalent  to  decreased 
bias)  when  adding  another  receiver  into  the  inverse  problem.  Resolution  index  curves  (proportional 
to  inverse  bias)  over  angle  are  compared  for  each  single-receiver  inversion  and  the  combined- 
receiver  inversion  case.  Note  that  the  benefits  seen  in  the  individual  receiver  results  are  essentially 
combined  when  both  receivers  are  used  -  each  receiver  provides  for  a  high  resolution  (low  bias)  in  a 
particular  angular  range,  which  together  provides  for  high  resolution  (low  bias)  over  a  much  wider 

angular  range. 


IMPACT/APPLICATIONS 


A  nonlinear,  well  characterized  filter-based  inversion  method  and  algorithm  will  have  application  to 
environmental  estimation  and  target  tracking.  A  practical  method  way  to  compute  the  resolution  for  a 
nonlinear  inversion  will  have  an  impact  on  the  characterization  of  uncertainty  and  uniqueness  of 
environmental  estimates  required  for  acoustic  propagation. 

RELATED  PROJECTS 

Our  research  is  directly  related  to  other  programs  studying  effects  of  uncertainty  in  the  environment, 
measurements,  and  models  on  acoustic  propagation,  and  target  detection  and  characterization. 
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